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ABSTRACT 


Photoelasticity is an optical method of experimental stress analysis, which yields a 

whole field representation of the principal stress difference and its orientation With 

advanced digital image processing systems, it has now become possible to get fringe 

order at every pixel in the domain by the use of phase map generation by interactive 

approach and phase unwrapping. Though for a large variety of problems phase 

unwrapping can be applied to whole domain in one step providing only one seed point, 

certain class of problems having complex geometiy need the individual processing of 
different zones. ^ 

The present work discusses the use of tiling for phase unwrapping. The purpose 
of the present work is to handle the phase map efBciently for conect unwrapping. The 
vanous difficulties related to digital implementations are discussed and have been taken 
care off. Various practical problems having complex geometo- are analyzed using 
completely experimental phase map extracted by six-step phase shifting technique. 

A comprehensive software in Windows environment was developed for 
photoelastic analysis. The software supports fringe order calculation by six step phase 
shiftmg algorithm, isoclinic calculation by Brown & Sullivan, Ramesh & Mangal and 

Pattemon & Wang algorithms. Stress separation by conventional Shear Difference and 

Tesar s Modified techniques is also incorporated. 
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CHAPTER 1 


INTRODUCTION 


1.1 Introduction 

Photoelasticity deals with the determination of stress fields by the use of polarized light. 
This method is based on the effect of temporary birefiingence. Certain transparent, non- 
cry stalline materials like plexiglass, celluloid, homolite and epoxy axe optically 
isotropic when not loaded but become optically anisotropic and display characteristics 
similar to crystals when they are stressed. The birefiingence in the material is retained 
only during the application of loads and disappear when they are removed. 
Measurements of this birefiingence provide data, sufficient enough for quantitative 
determination of the corresponding state of stress. The measurements are made with the 
help of an optical instrument called polariscope. A polariscope can be used to record 
maximum shear stress i.e., the difference in the principle stresses and their orientation. 
Corresponding contours are called isochromatic and isoclinic fiinges respectively. The 
experimentation in photoelasticity involves aquisition of stressed model images at 
various positions of optical elements, depending on the requirement of a specific 

algorithm. The main merit of photoelasticity is its simplicity in experimentation as well 
as in evaluation of stress field. 

The main advantage of photoelasticity is that, it can give the stress distribution 

for the whole field. In conventional methods tedious compensation methods are needed 

to obtain the fi'actional Singe order at any point other than the fringe position. Number 

of data points that could be considered are also less. These procedures are quite 

involved and requires skill in identification of isochromatics. Hence, automation of Hatg 

q isition and analysis to minimize these problems and to provide fast and accurate 

suits has become essential. With advanced digital image processing system, 

onsiderable success has been achieved in automating the data aquisition firom entire 

field and it has now become possible to evaluate the total fiinge order at every pixel in 
the domain. 

Photoelasticity directly provides only the information of principle stress 
difference and the direction of the principle stress at the point of interest. Using this, one 
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can find the normal stress difference and in-plane shear stress by invoking equations in 
mechanics of solids or Mohr’s stress circle. Determination of components of stress 
tensor requires the use of auxiliary methods. 

1 .2 Present Work 

A complete photoelastic analysis in a broad sense can be divided into five major steps, 

1 . Image aquisition 

2. Calibration of model material 

3. Retardation calculation 

4. Evaluation of isoclinic parameters 

5. Stress separation 

The software developed so far in the laboratory have remained in text based 
interface "with DOS environment. One of the major thrust of this thesis has been to 
convert the time tested software to Windows environment. 

In digital photoelasticity calibration involves Singe skeletonization, Singe 
thinning, data collection and analysis of these data in a least squares sense for obtaining 
the material sSess Singe value. Retardation calculation involves two major steps 
namely, phasemap generation and phase unwraping. Newer initiatives for phase 
imwrapping of arbitrary model geometry has been attempted m this thesis. 
Consideration of arbiSary model geometry requires phase unwraping in steps and calls 
for a tiling concept. The algorithms reported by Brown & Sullivan [1] and Ramesh & 
Mangal [2] (both based on plane polariscope arrangement) are fotmd to give good 
results [3] for isoclinic determination and are implemented in windows environment in 
the present work. Finally sSess separation using Shear Difference and Tesar’s modified 
technique [4 - 6] are also implemented for the determination of Ox , CTy and Xxy. 

1.3 Thesis Layout 

Chapter 2 provides an overview to data aquisition techniques. First histogram 
equalization and fringe thinning for image enhancement are discussed followed by 
boundary identification and its representation. Then various algorithms for fringe order 
and isoclinic determination are briefly discussed. 

Parameters affecting phase unwrapping [3] and use of tiling concept for phase 
unwrapping is discussed in detail in chapter 3. Phase unwrapping for the problems of 
disc and ring under diametral compression, rail road contact and spanner tightening the 
nut is discussed in detail. 
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Chapter 4 discusses the various problems in the evaluation of isoclinic 
parameter [7], The use of median filtering for improving the isoclinic result is 
discussed. Shear difference, Tesar’s stress separation techniques and their results for the 
problem of disk under diametral compression are discussed. 

Conclusion and scope for future work is presented in chapter 5. 



CHAPTER 2 


OVERVIEW OF DATA AQUISITION TECHNIQUES 


2.1 Introduction 

Isochromatics and isoclinics are the two fringe contours obtained in a photoelastic 
experiment. Isochromatic are contours of constant principle stress difference and 
isoclinics are the loci of points along which the principle stress direction (with respect to 
the reference axis) is a constant. For stress analysis, the photoelastic data required at any 
point are, 

• The isochromatic fringe order and 

• The isoclinic parameter at that point. 

The phase shifting technique for photoelastic analysis use direct processing of 
intensity data for quantitative extraction of isochromatic and isoclinic parameters for 
every pixel in the domain. Although in certain class of problems it may still be 
convenient to process a dark or bright field image to collect data of isochromatic fiinge 
order for a few selected points. Generally, the Singes do not appear as thin lines, but 
they appear as broad bands. So the first step in fringe processing is to identify the actual 
fringe in this dark band, i.e., fringe skeletonisation. 

2.2 Histogram Equalization [3] 

Histogram equalisation alters the dynamic range and contrast of an image. To perform 
histogram equalisation, one has to select a function representing a curve to modify the 
histogram. The available curves include Linear, Bell, Cube, Log and Exponential. This 
operation has the effect of lightening or darkening an overall image, depending on the 
curve chosen. The Linear equalisation curve moves the indexes so that the average 
index density across the grey-scale is even. This equalisation produces high contrast 
image with the highest possible dynamic range. The Bell equalisation curve moves the 
indexes in a similar manner, but groups the indexes more closely to the centre of the 
scale. This equalisation produces a high contrast image but not quite as high a dynamic 
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range as would achieved with a linear equalisation curve. The Cube and Log 
equalisation curves concentrate the range of indexes at tde low end of the scale, making 
the overall image darker. These equalisations produce a high contrast image "with a low 
dynamic range. On the other hand the exponential equalisation curve concentrates the 
range of indexes at the high end of the scale, making the overall image lighter. The 
modified image has a high contrast with low dynamic range. 

2.3 Thresholding [3] 

Thresholding is one of the most important technique to image segmentation and is 
essentially a process of separating background and object modes. The approach is to 
select a threshold T that separates these modes. Then, any point (x,3^) for which g{x,y)>T 
is called an object point; otherwise the point is called a background point. In a generic 
sense thresholding may be viewed as an operation that involves test against a function T 
of the form, 

T = f{x,y,p{x,y),g{x,y)} (2.1) 


Where g{x,y) is the gray level of point {x,y) and p{x,y) denotes some local 
property of these point, for example, the average gray level of a neighborhood centered 
on {x,y). A thresholded image b(x,y) is defined as. 


b{x,y) 


1 if gix,y)>T 
0 if g{x,y)<T 


( 2 . 2 ) 


When T depends only on g(x,>'), the threshold is called global. If T depends on 
both g{x,y) and piA,y), the threshold is called local. If, in addition, T depends on the 
spatial co-ordinates x and>', the threshold is called dynamic. 

2.4 Fringe Thinning 

Due to the nature of formation of fringes, the intensity of pixels in the fringe band 
(dark/bright) does not remain constant but a variation does exist. Though the variation is 
not resolvable to a fine degree by a human eye, the present day CCD cameras can easily 
recognize and quantize it for further processing. 

Although the use of intensity information leads to a better approach for fringe 
thinning, the early methods were confined to treating the fringe field as a binary image. 
Later, intensity based methods were developed. These techniques can be further 
classified into mask based methods and global methods. Since the global method is easy 
to implement and give better results, this was implemented in the software. 
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2.4.1 Global Identification of Fringe Skeleton Based on intensiti^ 
Variation 

In general, photoelastic images have high contraist and the edges of the firinges- caui Ibe 
extracted easily. Further the use of edge informal on to identify the fiinige skeleloffis 
eliminates the occurrence of spurious fringes. Tlray in global methods, fringe thirmiag 
involves two steps of edge detection and fringe skeJirtonization. 

Edge detection is an important step. The biestter is the edge detection algiori thm; 
the better will be the skeleton identification [8]. Eoige detection by thresholding iwoilks 
very well for photoelastic images. The algorithm is simple and also computationariy 
very fast. If there is intensity variations in the iimage, one can use either iynamic 
thresholding or divide the image into convenient ssegments and use appropriate gLoh al 
threshold for each of the segment to detect fringe ledges. The image thus obtained is tto 
be stored in one of the available frame buffers. 

The second part of the processing sequence is the skeletonization of the fringess. 
For this, fringe areas axe to be identified first and wdthin the fringe area, the s keLeton 
point needs to be scanned appropriately. One of simiplest scanning schemes is xo*AV-'wise 
scanning of the image and is applicable for fringes that are vertical or almost verticaaJ. 
The steps involved are, 

1. Scan the thresholded image row-wise from tlhe first pixel to identify the fringe 
edge. 

2. The scanning is carried out from one edge oftlke fringe to the other edge. Store tbie 
gray levels of the pixels within the boundaries ‘of the fringe. 

3. Identify the minimum intensity point and mark: is as a fringe skeleton pixel. 

4. Step (1) and (2) are repeated until all the pixefcs in the image are scanned and at tie 
end of the process, the thinned image is obtainaed. 

On the other hand, if the scanning is made? column-wise, the algorithm cam be 
applied to extract fringe skeletons which are prim«rily horizontal. To process coimptex 
fringe patterns, one requires a more refined scanoimg method taking into acco unt tie 
curvature and normal of the fringes. Though the i«Sea is conceptually very simple, it: is 
very difficult to implement from software point of view and is also computationallly 
very expensive. 

Fringe skeletons that are continuous and ftete of noise, can be easily extra.cte d by 
appropriately processing the results of fringe skelotons obtained by a finite number oof 
orthigonal scans. Once the edges are determined, the image is scanned row-wise (□“ 
scan), diagonal-wise (45° scan), column-wise (90° scan) and cross diagonal-wis e ( 135“ 
scan). Thus at the end of the initial processing, one^gets four images of fringe slceletoins 
corresponding to each scan direction. Each of these fringe skeleton is not contiiruoujs 
because of its direction dependence. The use of lotgical operators as shown in Fig . 2 -1 
can eliminate this discontinuities and noise [9~10], The logical “OR” operation behweea 
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the orthogonally scanned images helps to get a continuous fringe but the image 
contains noise. The noise is scan direction dependent and is removed by the logical 
“AND” operation performed between the logically “OR”ed images. 



Fig, 2.1. The scheme of logical operators to obtain continuous fringe-skeleton 


2.5 Boundary Identification 

To limit the determination of field parameter within the region of interest (within the 
specimen), the image has to be segmented. Segmentation of the image is done by 
detecting the boundary. In DIP literature, boundaiy of the image is detected using 
spatial filters like Robert’s, Prewitt, Sobel etc. These filters detect the edges based on 
abrupt change in gray levels. The edge detection algorithms not only detect the 
boundary, but also detect the regions in the image where there is abrupt change in gray 
levels. The set of pixels obtained from the above filters seldom characterises a boimdary 
completely. This is because, non-uniform illumination causes breaks in the boundary. In 
photoelasticity the accurate detection of the boundary is critical to apply algorithms for 
stress separation. Hence, instead of using the above general edge detection algorithms, 
the geometric properties of model can be used to detect the boundary accurately. In 
bright field arrangement of the polariscope, the boundary of the image can be seen 
clearly. Using this image, the boundary can be detected based on geometric primitives 
such as line, circle, ellipse and arc etc. The boimdary is to be drawn using primitives 
developed, such that outer side pixels of the line/curve drawn represents the boundary of 
the specimen. For ease of book keeping it is better to assign a specific gray level to the 
detected boundary. Once the boundary is identified, to make further computations, the 
boundary data has to be expressed in a suitable representation. The positions of the 









8 


pixels, which have the intensity of 255, are written into two separate intermediate files, 
one for horizontal scanning {scanx.dat) and one for vertical scanning (scany.daf). In 
view of spatial quantization, to represent the segment of the boundary having a 
curvature, more than one pixel is needed. Further processing of these files are done in 
which the first and the last column for every line for these files are extracted and stored 
into two separate files for horizontal scanning (*.ybn) and for vertical scanning (*jcbn). 

2.6 Plane Polariscope 

Plane polariscope is one of the simplest optical arrangement possible. Fig. 2.2 shows the 
optical arrangement for plane polariscope. A stressed model is kept in the field of view. 
The incident light on the model is plane polarised. As it passes through the model, the 
state of polarization changes from point to point depending on the principle stress 
direction and the magnitude of principle stress difference. The information about the 
stress field can be obtained if the state of polarization of the emergent light is studied. 
This is easily achieved by introducing a analyzer at 0°. With the introduction of 
analyzer, fringe contours that are black appear on the screen. The fringe contours 
correspond to those points where intensity of light transmitted is zero. Since the 
analyzer is kept at 0°, this is possible only if the emergent light from the model has its 
plane of polarization vertical or in other words the incident plane polarized light on the 
model is unaltered. 



Fig. 2.2 Generic arrangement of a plane polariscope 

The incident light is unaltered on all those points where the model behaves like a 
full wave plate. This happens when the principle stress difference (ari-cy 2 ) is such as to 
cause a relative phase difference of 2m7t (m = 0, 1, 2, ...), where m is an integer. Since 
stress is continuous one observes a collection of points forming contours that satisfy the 
above condition and the respective fringe field is known as isochromatics. 
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Another possibility, wherein the incident light is unaltered is, when the 
polarizer axis coincides with one of the principle stress directions at the point of 
interest. In this case, light extinction is not wavelength dependent and one observes a 
dark fringe even in white light. These are known as isoclinics meaning contours of 
constant inclination. The principle stress direction on all points lying on a isoclinic is a 
constant. 

2.6.1 Intensity of Light Transmitted for a Generic Arrangement of a Plane 
Polariscope [3] 

Figure 2.2 shows a photoelastic specimen kept in a plane polariscope with the 
polarizer and the analyser kept at arbitrary angles of a and /? respectively. Using 
Jones calculus, the components of light vector along the analyser axis and 
perpendicular to the analyser axis are obtained as, 


f 
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The intensity of light transmitted /p (the subscript ‘p’ denotes that the incident 
light is plane polarised) is obtained as. 




5 5 

cos^— cos^(y5 -a) + sin^— cos^(y9 + a- 29) 


(2.4) 


where 4 accounts for the amplitude of incident light vector and the 
proportionality constant. 


2.7 Circular Polariscope 


In the case of a circular polariscope, a circularly polarized light is used to reveal the 
stress information in the model. This is achieved by introducing a quarter wave plate 
with one of its axis at 45° after the polarizer. The emergent light is analyzed using a 
combination of quarter wave plate and an analyzer. The basic arrangement is shown in 
Fig. 2.3. Inside the model boundary, the intensity of light transmitted is governed by the 
stress field and this can be obtained by Jones calculus. 



2.7.1 Intensity of Light Transmitted for a Generic Arrangement of a 
Circular Polariscope [3] 
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Fig. 2.3. Generic arrangement of a circular polariscope 

Consider the optical arrangement shown in Fig. 2.3 in which the second quarter-wave 
plate and the analyzer are kept at arbitrary positions. The first quarter-wave plate is 
kept at = 135-deg. Using Jones calculus, the components of light vector along the 
analyzer axis and perpendicular to the analyzer axis (for 135-deg) are obtained as, 


E, 1 1 

COS P 

sin p 

[■^jS+r/2 J 2 

-sin p 

cos P 


l-icos 27 -isin 27 
i sin 2 ^ 1 + i cos 2t] 

s s s ~ 

cos isin— cos 2 ^ -isin— sin 2 ^ 

2 2 2 

5 S S 

-isin— sin 2 ^ cos— + i sin— cos 2 ^ 

2 2 2 . 



\m 


(2.5) 


The combination of the quarter- wave plate and the polarizer for ^ = 135-deg 
produces a left circularly polarised light. Denoting the intensity of light transmitted as I^, 
the intensity expression is obtained as, 

/j = ^ ^ [sin 2{p - rj) cos 5 - sin 2{6 - rf) cos 2 (/? - 77 ) sin <J] ( 2 . 6 ) 


where Ja accounts for the amplitude of light vector and proportionality constant. The 
intensity equation developed show that when the extinction conditions are satisfied, the 
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intensity of light transmitted should be zero. However, this never happens in an 
experiment. In order to account for this behaviour a term /b is added to account for the 
stray light/background illumination. 

In subsequent discussions, the intensity of light transmitted will be denoted as If 
for the optical arrangement. The orientation of the optical element that needs to be 
rotated will also appear with a subscript i. Thus, in Eq. (2.6) by adding the term 4 to 
account for background illumination/stray light reduce to 

/; ~ ^ - 77 , ) cos 6 - sin 2{9 - 77 ,) cos 2(/?, - 7 j) sin J] (2.7) 

2.8 Six-step Phase Shifting Algorithm [11] 


One of the most successful and time-tested approach of phase shifting to evaluate 
isoclinic and isochromatic parameter involves the use of six different optical 
arrangements. The various positions and the respective intensity equations are 
summarized in Table 2.1. The arrangement shown also accounts for the mismatch of 
quarter wave plate error. 

Using the intensity equation listed in Table 2.1, the isoclinic parameter 9 is 
obtained as. 
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( 2 . 8 ) 


From the equations shown in Table 2.1, the fractional retardation 5is obtained as. 
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(2.9) 

( 2 . 10 ) 


In Eqs. (2.8) to (2.10) the subscript ‘c ’ for the 0 or 0 indicates that the principal 
values of the inverse trigonometric function is referred. In subsequent discussion too the 
same symbolism will be used. 

It is clear that the exact calculation of isoclinic angle and fi-actional retardation is 
not trivial from Eqs. (2.8) to (2.10). Equation (2.8) shows that the isoclinic value is 
indeterminate at those points where 0(, “ 27 r, Stt .... The implication is that the 

isoclinic contour will not be continuous over the domain. When 5^ is not exactly equal 
to 0 or TT or 2n etc. but very close to these values, 9^ does not become indeterminate but 
the value of 0 c determined will be unreliable and appears as noise in isoclinic fringe 
when plotted. This is known as isochromatic-isoclinic interaction in the calculation of 
the isoclinic parameter. 
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Table 2.1. Polariscope arrangement and intensity equations for six-step phase shifting method. 
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Intensity equation 
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1 3 m 

1 

m 

\ 

4 = 4 + ^ (1 “ cos 29 sin S) 


Considering that 6^ is calculated correctly by Eq. (2.8), the calculation of by 
Eqs. (2.9) or (2.10) should result in the same value of if continuous variables are used 
in the arctangent calculations. In practical applications, this is not the case due to 
intensity digitisation and quantization and only discrete values are used in the arctangent 
calculations. This can lead to error in the calculations. Further, in Eqs. (2.9) and (2.10), 
has low-modulation regions when cos2^c or sin20c are small. In these regions, the 
intensity differences (I4 - 4) or (/; - 1^) will also be small. The evaluation of will 
involve the ratio of two small discrete numbers and will be in considerable error in such 
regions. This is termed as isoclinic-isochromatic interaction in the calculation of the 
fractional retardation. Apart from this, the principal value of evaluated by Eq. (2.8) 
lies in the range -jt/4 to 7r/4 and an ambiguity exists on whether 9^ corresponds to cTj or 
CTj direction over the domain. 

It is instructive to note that regions of low values of cos2^c coincide with high 
values of sin2^c and vice versa. Therefore for every point in the region of interest one 
has two estimates of 5^ with different degrees of error depending on the region. Thus, 
the results must be combined appropriately to get the best possible values of 5^. For the 
whole domain it is obtained as. 


= tan" 
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= tan~* 

f 4sin^- 
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( 2 . 11 ) 
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(h-h) 

Unlike Eqs. (2.9) or (2.10), Eq. (2.1 1) always guarantees a high modulation over 
the field [12]. 
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Figure 2.4 shows the sequence of six images recorded for a disk under 
diametral compression (dia 60.10 mm, thickness 6.28 mm and load 773.0 N). The 
material stress fringe value is calculated by digital technique and is found to be 1 1.9973 
N/mm/fr. Figure 2.5a shows the phase map obtained by Eq. (2.1 1). on scrutiny one can 
see that near the load application points the phase map is not good. This is not just due 
to high stress gradient but is due to the interaction of isoclinic on the isochromatic 
evaluation. The principal values of 0^ of Eq. (2.8) lies in the region -rt/A < 0^ < ^/4 
whereas, physically, principal stress direction lies in the range -7tl2 < 0 < 7tj2 and this 
introduces an ambiguity on 5 evaluation. Figure 2.5b shows the phase map obtained 
from experimental images using the theoretically calculated Vindicating the cTj direction 
throughout the domain. The phase map is uniformly good over the domain. 

Though, phase shifting can provide 5 over the whole domain, certain level of 
human interaction is needed to correctly identify its sign for unwrapping. This is one of 
the serious limitation of the technique and could be resolved only when isoclinic is 
unwrapped correctly over the domain corresponding to either cr, or <72 direction 
uniformly over the domain. This problem is current area of research and automatic 
solution is still elusive. Figure 2.5c shows the phase map of 2.5a improved by the 
interactive approach developed by Ramesh and Mangal [13]. The approach is general 
one and is applicable for any problem and one requires the information of fringe 
gradient in the bad zone for phase map improvement. Figure 2.5d shows another way of 
phase map improvement which is specific to the problem of circular disc under 
diametral compression by unwrapping the Vc in the following maimer. In the first and 
third quadrant , 


0, if(V>0) 

0^+nl2 if(V<0) 


and in the second and fourth quadrant , 


Vc if (V < 0) 

9^- It 11 if(V>0) 


( 2 . 12 ) 


(2.13) 


The basic advantage of phase shifting is that inverse tangent functions are used 
for 5 calculation. The unvwapping of <5 becomes a lot more simpler and straight forward. 
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Fig. 2.4. Sequence of six images recorded of a disk under diametral compression (dia 60.10 mm, 
thickness 6.28 mm and load 773.0 N) for the optical p>osition listed in Table 2.1 
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Fig. 2.5. Whole field representatioa of fractional retardation, 4 in the form of phase map of a disk under 
diametral compression obtained using Eq.(2.il) with 4 being calculated by various means: a 4 by Eq. 
(2.8) and 4 by Eq. (2.11) b 4 by the*ory corresponding to the cr^ direction and 4 by Eq. (2.1 1) c phase 
map Fig. 2.3a improved by interactive approach d 4 by Eq. (2.8) then unwrapped by Eqs. (2.12) & (2.13) 
according to quadrant and 4 by Eq. (2.11) 


2.9 Evaluation of isoclinic parameters 

There are a number of algorithms available for isoclinic angle calculation based on 
either circular or plane polariscope arrangement. The algorithm in plane polariscope 
arrangement are, 

• Ramesh and Mangal [2] 

• Brown and Sullivan [1] 

• Chen and Lin [14] 

• Dupre et al. [15] 

• Sarmaetal. [16] 

• Kihara [17] 



16 


Of these, the methods of Brown & sullivan and Chen & Lin are actually 
called as polarization stepping methods while the remaining methods are called as phase 
shifting methods. 

Out of these algorithms, the algorithm given by Ramesh & Mangal and Brown 
& Sullivan are found to give good results [7] and are implemented in windows 
environment. 

2.9.1 Aigorithm by Ramesh and Mangal [2] 

In this, the polarizer and analyzer are kept at arbitrary positions. The intensity of 
light transmitted is given by Eq. (2.4). If the term 4 is added to account for stray 
light/background illumination, one gets 




cos^— + sin^— sin^( 20 -/?, ) 


(2.14) 


If one takes all possible combinations of cc and /? in steps of ^/4, then only six 
intensity equations turn out to be unique which are given in Table 2.2. Out of these 
only the first four intensity equations, the isoclinic parameter is calculated as. 
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Use of atan2 function is recommended for evaluating 0^ Eq. (2.15). Figure 
2.6 shows the experimentally recorded images for the disk. 

2.9.2 Algorithm by Brown and Sullivan [1] 

In this polarizer and analyzer are kept crossed and four images are obtained at step of 
7 r /8 radians. The polarisation steps used by Brown and Sullivan and the appropriate 
intensity equations are summarised in Table 2. 3. Figure 2.7 shows the experimentally 
recorded images for the ring. Isoclinic parameter is calculated as. 
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Use of atan2() ftmction is recommended for evaluating above equation. 
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Tabie 2.2. Table showing the intensity equations for different orientations of the polarizer and 

analyzer 

Ramesh and Manga! Algorithm 
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Table 2.3. Polariscope arrangements and intensity equations for polarization stepping schemes. 
Brown and Sullivan (Polarizer & Analyzer crossed; a=/?+;r/2) 
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Fig. 2.6. Sequence of six images recorded of a disk under diametral compression for the optical position 
listed in Table 2.2 
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2.10 Closure 

An overview to data aquisition techniques was presented in this chapter. Histogram 
equalization and fringe thinning algorithm are presented. Six step phase shifting 
algorithm for isochromatic determination was discussed. Algorithms reported by 
Ramesh & Mangal and Brown and Sullivan for isoclinic determination were discussed 
in brief 


Fig. 2.7. Sequence of four images recorded of a disk under diametral compression for the optical position 

listed in Table 2.3 



CHAPTER 3 


PHASE UNWRAPING BY TILING 


3.1 Introduction 

In the previous chapter, six-step phase shifting algorithm to obtain fractional retardation 
over the complete model domain has been presented. In this technique, the fractional 
retardation is represented as a phase map. In the phase map, the fractional fringe orders 
in the range of 0-1 are represented as an intensity map of 0-255 gray level. For practical 
utilization of the data, one has to find the total fringe order over the domain. This is 
achieved by a process called phase unwrapping. For phase unwrapping to be effective, 
phase map should be free of noise and discontinuities. 

Though the algorithm for phase unwrapping is simple and straight forward, there 
are certain practical problems that affect phase unwrapping. Selection of threshold gray 
level Tp, for un\vrappmg and location of primary seed point are two major parameters 
that affect the phase unwrapping process [3]. 

Problems emerging from choice of primary seed point can be overcome by the 
use of tiling concept, and is discussed in detail in this chapter. With the use of tiling, the 
area of phase map containing noise can be treated separately and propagation of error 
due to noise can be avoided. 

3.2 Noise Removal in Phase Maps 

The phase maps in general may have salt and pepper noise (binary noise) at a few 
points in the domain. To remove noise from images which have sharp edges, usually a 
median filter is used as other filters can destroy the sharp edges. The median filter is a 
non-linear spatial filter. In this, the grey level of each pixel is replaced with the median 
of the grey levels in the neighbourhood of that pixel. Even though, the median filter 
preserves the sharp edges, some information at the sharp edges is lost. The intensity 
gradient at the transition of wrapped regions is the main parameter for unwrapping of 
the phase. Any decrease in the intensity at the transition regions will affect phase 
unwrapping adversely. Thus, median filter is not an ideal choice for the present 
application. For removing the noise without disturbing the sharp edges, the following 
noise removal algorithm can be used. 
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This algorithm removes the binary noise based on the difference in the 
intensities between the neighbourhood pixels along the scanning direction. If the 
difference in intensities between the current pixel and the neighbourhood pixels along 
the scanning direction is beyond a threshold, then the current pixel is identified as a 
noise pixel. The grey level of this pixel is replaced with the average of grey levels of the 
neighbour-hood pixels. 

Let the noise removal gradient threshold be T„. Scanning the image horizontally, 
the pixel at (0, 0) is a salt noise pixel when the following condition is satisfied. 

(Soo - g-!o)> T„ and (gio - goo)< -Tn (3.1) 

The pixel at g(0,0) has pepper noise w'hen the following condition is satisfied. 

(goo-g-!o)> -T„ and (gio-goo)<T„ (3.2) 

The grey level of the noise pixel g(0, 0) is then replaced with the grey levels of 
the neighbourhood pixels [(g-io + gio)/2]. In practice, the noise removal scheme could be 
merged with the program doing fractional retardation calculation itself. 

The value of r„ depends on the maximum grey level gradient between the 
neighbourhood pixels in the data region and one has to select it by trial and error. 

3.3 Algorithm for Phase Unwrapping [3] 

Figure 3.1 shows the phase map for the case of a disk under diametral compression. The 
variation of intensity in the phase map along the horizontal diameter is shown in Fig. 
3.2. It can be seen that the fractional fringe order increases in the direction in which 
total fringe order increases and a discontinuity in the slope occurs when the phase 
retardation \s 2 tt. 



Fig. 3.1. Typical phase map for the problem of disk under diametral compression 
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xtR 

Fig. 3.2. Fringe order variation along the horizontal diameter for the problem of disk under diametral 
compression 

If an integer fringe order is added or subtracted at the transition of the wrapped 
regions, the total fringe order can be known. The slope of the curve can decide addition 
or subtraction. If it is positive, then, one has to add the integer value, otherwise subtract. 
The methodology however requires the total fringe order for at least one point in the 
domain to be supplied by atixiliary means. The point at which total fringe order is 
interactively supplied is termed as primary seed point. 

To unwrap the phase completely within the boundary, the image has to be 
scanned completely. A simple procedure is to scan the image along the row passing 
through the primary seed point followed by column-wise scanning of the entire image. 
The initial scanning along the row provides seed points for effecting the scanning along 
the columns. These seed points are termed as secondary seed points. The entire scheme 
can be termed as row- wise seeding and column-wise scanning. 

Let point ‘A’ in Fig. 3.2 be the primary seed point at which the total fringe order 
is to be supplied externally. Once the point is selected the software asks for unwrapping 
threshold Tp and total fringe order at that point. The unwrapping is first done along the 
left direction and then along the right direction. At point ‘B’ there is a discontinuity in 
the slope. This discontinuity in the slope is identified by the threshold value Tp . When 
the slope between the neighborhood pixels is less than -Tp , the total fringe order is 
increased by one and when the slope between neighborhood pixels is greater than Tp , 
the total fringe order is decreased by one. 

Though in principle the fractional fringe order gradually increases in the 
direction of the overall fringe gradient, there can be kinks in the curve, due to random 
noise, which escaped the noise removal procedure. Since phase unwrapping is basically 
a scaiming procedure, this noise can easily propagate and will adversely affect the final 
results unless corrective measures are taken. Every effort must be made to select the 
primary seed point such that, no noise is present in the row containing pr imar y seed 
point. 
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The other scanning algorithm first generates secondary seed points in the 
column containing primary seed point followed by row-wise scanning. This is called 
column-wise seeding and row-wise scanning. 

3.4 Parameters Affecting Phase Unwrapping 

3.4.1 Influence of the Selection of the Phase Unwrapping Threshold [3] 

Unwrapping threshold identifies the transition regions of phase map. In analog data 
processing one can take as 255. However, in digital data processing even in small 
stress gradient zones, one has to take a smaller value of than 255 to accommodate the 
effect of spatial discretization and intensity quantization. This aspect is prominent in 
high stress gradient zones. Variation of the intensity in phase map along vertical 
diameter for the problem of circular disc under diametral compression is shown in Fig. 
3.3. It can be observed, that, at stress concentration regions, the intensity in the wrapped 
regions neither starts with a value of 0 nor ends with a value of 255. This is primarily 
due to the fact that the fringe gradient is so high that one requires sub-pixel calculations 
to identify the exact transition zone. 



Fig. 3.3. Variation of intensity in phase map of Fig. 3.1 along the vertical diameter 

The error in mismatch of quarter-wave plates also plays a role in deciding 
the value of T^. Error in mismatch of quarter-wave plates has the effect of slightly 
shifting the point of transition and the intensity does not reach the maximum value 
of255 [18]. 

From the above discussion it is clear that should be lower than 255. In 
general, should be higher than the threshold, that is used to identify noise in 
the domain. 

3.4.2 Influence of the Location of the Primary Seed Point [3] 

To initiate the phase unwrapping procedure, the user has to specify the total fringe order 
for at least one point in the field. Ideally, the location of this point should be arbitrary 
and the phase unwrapping procedure should be flexible enough to unwrap the phase for 
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the entire domain. However, the phase unwrapping algorithm could be a lot simpler if 
the user selects the primary seed point appropriately. The choice of location of seed 
point depends largely on the specific geometry of the specimen under consideration. 

For example, for the problem of disk under diametral compression, while 
using the row-wise seeding and column-wise scaiming scheme to unwrap the entire 
image, the primary seed point should be chosen along the horizontal diameter. If 
not, the secondary seed points will be generated only up to the boundary, that is 
encountered along the row-wise seeding. On the other hand, if column-wise 
seeding and row-wise scanning is used, the primary seed point must lie on the 
vertical diameter. In general, the primary seed point should be located on a 
horizontal line covering the entire width of the model, if row-wise seeding and 
column wise scanning is used and vice versa. 

This problem is further enhanced if the model geometry is complex. For 
example, for the problem of annular ring under diametral compression, where no 
single primary seed point can be specified which will efficiently unwrap the phase 
map. In this case, if true boundary information is used then there will always be 
limited seeding. Otherwise, if only outer boundary information is used, for all 
points inside the inner boundary, the secondary seed points generated will not have 
any significant secondary seed value, as phase map itself contains noise in that 
zone. This problem of limited unwrapping and insignificant seeding can be 
overcome by adopting a tiling procedure for phase unwrapping. 

3.5 Use of Tiling Concept for Phase Unwrapping 

Tiling in a broad sense can be understood as, treating the whole image as an 
assembly of small windows called tiles. The process parameters can individually 
be selected for each of the tiles and processed information of each tile can then be 
reassembled to obtain the final results (or overall image) as if the whole image was 
processed simultaneously. In image processing, such operations are mostly used 
for intensity manipulations and can be easily accomplished by using the frame 
buffer memory. However, in case of phase unwrapping this approach is not 
possible as unwrapped fringe orders are real numbers, that can’t be stored in pixel 
locations of a frame buffer. For storing the fringe order values, a file structure has 
to be used. This poses certain restrictions on the selection of the nature of the tiles. 
Selecting an arbitrary window would be difficult from book keeping point of view. 
Even if an efficient file structure is used, pasting of processed information back to 
the original place from each of the file is a difficult process and will require 
tedious programming efforts. This restricts the shape of tiles either as vertical 
strips filling the entire height of the image or as horizontal strips filling entire 
width of the image. 
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An alternative approach is to declare a large matrix of real numbers, 
which is equal to size of the image. In such case any arbitrary window can be 
selected. Though this is not a feasible approach on the systems having small RAM 
size, but for present day systems with RAM size exceeding that of 64 Mb, this 
approach is quite useful and efficient. Also maintaining a file structure, as required 
in first approach is logically difficult, slow in operation and more error prone. In 
the second approach, since the data are dynamically read from RAM and written to 
RAM, operations are quite fast. 

Figure 3.4 shows typical problem of a plate with a hole, subjected to uniaxial 
tension. Due to the symmetry of the problem one quarter of the plate could be 
considered for evaluation of the total fringe order. If the row-wise seeding and column- 
wise scanning is used, fringe order in the complete domain will be successfiilly 
xmwrapped only if the primary seed point lies in the zone 1 of the Fig 3.4a. Conversely 
if column-wise seeding and row-wise scanning is used, primary seed point should lie in 
the zone 1 of Fig. 3.4b. If this is done phase map will be properly unwrapped. 

If for some practical reasons, it is not possible to locate the primary seed point in 
these zones, viz. seed point is selected in zone 1 of Fig. 3.4a and row-wise seeding is 
used. Since in zone 1, width of model is smaller than in zone 2, this would lead to 
limited seeding only up to left boundary and secondary seed value for vertical scanning 
will be unreliable outside the boundary. Similarly, if primary seed point in Fig. 3.4b is 
selected in zone 2 with column-wise seeding and row-wise scanning then secondary 
seeding will be limited only up to lower boundary of zone 2 and problem of limited 
seeding arises again. 

In such cases use of tiling could help in imwrapping the complete domain 
effectively. For row-wise seeding and colunrn wise scanning if the primary seed point is 
in zone 2, then, the tiles needed for unwrapping are shown in Fig. 3.4c. Tile 1 will be 
unwrapped first based on primary seed value specified at point ‘A’, followed by row- 
wise seeding and column-wise scanning in zone 2, with primary seed value specified at 
any point on jimcture of tiles. Since tile 1 is already unwrapped seed point for tile 2 can 
easily be located at any point on the common botmdary of tile 1 and tile 2. Similarly for 
column wise seeding and row- wise scanning the tiles needed are shown in Fig. 3.4d. 

Besides this, the use of tiling avoids error propagation due to noise, e.g. if phase 
map contains noise in some zones and scanning crosses these zones. The error due to 
noise of these zones will propagate and unwrapped fringe orders become unreliable. On 
the other hand, if tiling is used, then the tiles can be arranged so that in no case the 
scanning is continued after noisy zones. In such cases, the good zone of phase map will 
be unwrapped separately as another tile and seed value for this tile will be 
mdependently supplied. 

A software has been developed implementing these ideas and it supports the 
use of separate unwrapping threshold value for different tiles. User can 
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interactively supply threshold value, along with the seed point in each of the 
tiles. The softw'are supports both, row- wise seeding with column- wise scanning and 
column-wise seeding with row-wise scanning. Any of these two can be interactively 
used depending on the location of particular tile. The user has to make the tiling plan 
through menu selection and the tiles selected are internally stored in a tiling information 
file named tile.dat. 



Fig. 3.4. Tiling for phase unwrapping. Identification of zones for locating the primary seed point: a row- 
wise seeding and column-wise scanning b column-wise seeding and row-wise scanning. Tiles required 
for c row-wise seeding and column-wise scarming d column-wise seeding and row-wise scanning. 

Some general guidelines for proper selection of tiles are discussed in subsequent 
topics through the problems of disk and ring under diametral compression, rail road 
contact and spanner tightening the nut. 

3.6 Tiles for Problem of Disk 

Due to symmetry of model the discussion is restricted to one quarter of the image. As in 
this problem supplying the primary seed point along the horizontal or vertical diameter 
is not a problem, so the use of tiling does not serve any major purpose. If due to some 
problems, viz. presence of noise in phase map along these diameters, it is not possible to 
select a seed point along these diameters then any of the tiling arrangements shown in 
2 nd, ^rd ^th q^^drant (Fig. 3.5a) can be used to effectively unwrap the phase map. 

In practice, the tiles chosen can overlap. In order to test the implementation of 
the software, tiling plans used for the problem of disk imder diametral compression are 
shown in Fig. 3.5b, 3.5c and 3.5d. Tables 3.2-3.4 give the corresponding unwrapping 
parameters. The phase map is shown in Fig. 3.5e. Figure 3.5f shows the location of 
points for which sample data of unwrapped fringe orders are collected for comparison. 
These points correspond to points noted in Table 3.5. 
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Fig. 3.5. a Some tiling arrangements for the problem of disk under diametral compression b,c,d Tiles 
used for umwapping corresponding to Table 3,2“3.4 e Phase map used for unwrapping f Sample data 
points marked on the dark field image corresponding to Table 3,5. 
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Table 3-1 Unwrapping peirameters used for unwrapping complete image for problem of disk 


Tile 

Tile Coordinates 

Seed Point 

Seed Value 
(N) 

Threshold 

(tp) 

Seeding 

1 

(0,0) (511,511) 

(256,170) 

2.00 

140 

Row-wise 


Table 3.2 Unwrapping parameters and tiling plan 1 for problem of disk 


Tile 

Tile Coordinates 

Seed Point 

Seed Value 
(N) 

Threshold 

(tp) 

Seeding 

1 

(7,9) (442,267) 

(343,251) 

2.00 

140 

Row-wise 

2 

(6,251) (367,501) 

(170,262) 

2.00 

140 

Row-wise 

3 

(340,88) (506,504) 

(402,263) 

1.00 

140 

Row-wise 


Table 3.3 Unwrapping parameters and tiling plan 2 for problem of disk 


Tile 

Tile Coordinates 

Seed Point 

Seed Value 
(N'k 

Threshold 

(tp) 

Seeding 

1 

(13,13) (259,257) 

(170,246) 

2.00 

140 

Row-wise 

2 

(8,248) (504,385) 

(111,254) 

1.00 

140 

Row-wise 

3 



1.00 

140 

Row-wise 

4 

(251,7) (507,253) 

(305,143) 

3.00 

140 

Column-wise 


Table 3.4 Unwrapping parameters and tiling plan 3 for problem of disk 


Tile 

Tile Coordinates 

Seed Point 

Seed Value 
(N) 

Threshold 

(tp) 

Seeding 

1 

(12,122) (67,413) 

(17,253) 

0.00 

140 

Row-wise 

2 

(57,67) (157,483) 

(111,255) 

1.00 

140 

Row-wise 

3 

(152,11) (366,500) 

(170,257) 

2.00 

140 

Row-wise 

4 

(360,24) (506,492) 

(402,254) 

1.00 

140 

Row-wise 

5 

(96,26) (154,75) 

(151,71) 

1.00 

140 1 

Row-wise 


Three dimensional plots of the unwrapped firinge order for each of the 
arrangements of Table 3. 1-3.4 are shown in Fig. 3.6a, 3.6b, 3.6c and 3.6d respectively. 
It can be seen that fringe order plots are not sharp near the point of application of load at 
upper end. This is because in that zone phase map contains noise. Also, in second tiling 
arrangement unwrapping in that zone is very bad and gives negative fringe orders. This 
is because the column- wise seeding used for that tile crosses noise in phase map. Thus 
secondary seeding is completely disturbed due to noise and subsequent row-wise 
scanning looses its significance. This indicates, that the bad choice of seeding and 
scanning directions can adversely affect the results. The results summarized in Table 3.5 
indicate the success in implementing tiling concept for phase unwrapping. 
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3.7 Tiles for Problem of Annular Ring 

A few tiling arrangements for the problem of the annular ring under diametral 
compression is sho\\n in Fig. 3.7a. It can be noted that the tiling arrangement is so 
chosen, that either of the row-wise or column wise seeding in each of tile will serve the 
purpose, if the primary seed point is selected properly. Since for all practical purposes it 
is difficult to judge the total fringe order in tiles 3 and 5, the unwrapping procedure 
should be started from any one of the 1®\ 2"**, 4* or 6* tile and then be extended to 3'^'* 
and 5* tiles. The tiles can be made overlapping to each other to further simplify the 
selection of seed points in subsequent tiles. 

Tables 3.6 to 3.8 show the different tiling arrangements and unwrapping 
parameters used to unwrap the phase map of Fig. 3.7d for the problem of annular ring 
under diametral compression with 512 x 512 image size. Corresponding tiles are shown 
in Fig. 3.7b and 3.7c respectively. Table 3.9 shows the unwrapped fringe orders at some 
salient points for each of these tiling arrangements. These sample points are marked on 
the phase map in Fig. 3.7e. Figure 3.8a shows the three dimensional plot for fringe order 
data for ring, generated theoretically. Three dimensional plots of the unwrapped fringe 
order for each of the tiling arrangements of Tables 3.6-3. 8 are shown in Fig. 3.8b, 3.8c 
and 3.8d respectively. 


Table 3.6 Unwrapping parameters used for unwrapping complete image for problem of ring 


Tile 

Tile Coordinates 

Seed Point 

Seed Value 
(N) 

Threshold 

(tp) 

Seeding 

1 

(0,0) (511,511) 

(372,111) 

3.00 

140 

Row-wise 


Table 3.7 UnwTapping parameters and tiling plan 1 for problem of ring 


Tile 

Tile Coordinates 

Seed Point 

Seed Value 
(N) 

Threshold 

(tp) 

Seeding 

1 

(45,16) (470,139) 

(376,119) 

3.00 

140 

Row-wise 

2 

(45,375) (467,495) 

(403,387) 

2.00 

140 

Row-wise 

3 

(11,139) (136,375) 

(123,185) 

3.00 

140 

Row-wise 

4 

(378,139) (502,375) 

(412,357) 

2.00 

140 

Column-wise 

5 

1 

1 

1 

a 


3.00 

140 

Row-wise 

6 

(256,256) (378,375) 

(338,368) 

3.00 

140 

Row-wise 

7 

(136,256) (256,375) 

(178,371) 

3.00 

140 

Row-wise 

8 

(256,139) (378,256) 

(377,166) 

3.00 


Column-wise 











31 



Fig. 3.7. a Some tiling arrangements for the problem of ring under diametral compression b,c Tiles used 
for unwapping corresponding to Table 3.7-3. 8 d Phase map used for unwrapping e Sample data points 
marked on the phase map corresponding to Table 3.9. 
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Table 3.8 Unwrapping parameters and tiling plan 2 for problem of ring 


Tile 

Tile Coordinates 

Seed Point 

Seed Value 
(N) 

Threshold 

(tp) 

Seeding 

1 

(254,17) (470,138) 

(379,125) 



Row-wise 

2 

(376,132) (505,264) 

(411,210) 

2.00 

140 

Row-wdse 

3 

(254,138) (376,254) 

(341,144) 

3.00 

140 

Row-wise 

4 

(136,138) (254,264) 

(151,162) 

3.00 

140 

Column-wise 

5 

(136,264) (254,374) 

(175,369) 


140 

Row-wise 

6 

(254,264) (376,374) 

(336,369) 


140 

Row-wise 

7 

(40,15) (259,144) 

(134,116) 

3.00 

[iimQiii[ii 

Row-wdse 

8 

(35,141) (138,298) 

(95,175) 

2.00 


Row-wise 

9 

(10,155) (52,390) 

(47,258) 

1.00 


Row^-wise 

10 

(46,293) (155,477) 

(129,341) 

3.00 

140 

Row-wise 

11 

(148,370) (476,498) 

(372,396) 

3.00 

140 

Row-wdse 

12 

(372,260) (504,377) 

(413,315) 


140 j Row- wise 


Table 3.9 Fringe orders at some salient points unwrapped with different tiling arrangements shown in 
Tables 3.6-3. 8 


Coordinates 



Unwrapped Fringe Order (N) 


Without Tiling 


2.62 2.62 2.62 



















































Fig. 3.8 Three dimensional fringe order plots: a Theoretical data for ring b,c,d phase map of Fig. 
3.7d unwrapped with different tiling arrangements shown in Tables 3.6-3. 8 respectively 

Note that for points 8-12 unwrapped fringe order with complete image are not 
matching back to corresponding phase map and are shown bold. Since magnitude of 
error is almost one complete iBringe order, it can be because of some noise pixel causing 
propagation of error to these zones. These fringe orders are improved after use of tiling 
as shown in 2"“^ and 3'^'* tiling plan fringe order column. 
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3.8 Tiling Applied to Rail Road Contact Problem 

Tiling arrangement for the left half of the model is shown in Fig. 3.9a. Here, all tiles can 
be merged into a single large tile as long as it is possible to give primary seed point in 
the zone 2 followed by row-wise seeding and column-wise scanning. Otherwise one 
should go for tiling. It is interesting to note that, if 4* and 5* tiles are merged and seed 
point is given at top or bottom boundary and row-wise seeding is done, the concave 
boundary of the model will be encountered twice in some of the vertical scans. This will 
create problem as the phase information outside the model will be unreliable and error 
will propagate. So it is better to keep 4*^ and 5‘*’ tile as two separate tiles. 

Tables 3.10 to 3.12 show the different tiling arrangements and corresponding 
imwrapping parameters used to unwrap the phase map of Fig. 3.9d for the rail road 
contact problem with 512x512 image size. Corresponding tiles on the dark field image 
are shown in Fig. 3.9b and 3.9c. Table 3.13 shows the unwrapped fringe orders at some 
salient points for each of these tiling arrangements. These points are marked on 
corresponding dark field image and shown in Fig. 3.9e. Three dimensional plots of the 
unwrapped fringe order for each of the tiling arrangements of tables 3.10-3.12 are 
shown in Fig. 3.10a, 3.10b and 3.10c respectively. 


Table 3.10 Unwrapping parameters used for unwrapping complete image for rail road contact problem 


Tile 

Tile Coordinates 

Seed Point 

Seed Value 
(N) 

Threshold 

(tp) 

Seeding 

1. 

(0,0) (511,511) 

(111,231) 

1.00 

140 

Row-wise 


Table 3.11 Unwrapping parameters and tiling plan 1 for phase unwrapping of rail road contact problem 


Tile 

Tile Coordinates 

Seed Point 

Seed Value 
(N) 

Threshold 

(tp) 

Seeding 

1. 

(32,211) (483,329) 

(160,246) 

2.00 

140 

Row-wise 

2. 

(44,82) (491,222) 

(108,184) 


140 

Row-wise 

3. 

(259,324) (391,437) 

(294,396) 



Column-wise 

4. 

(118,326) (263,441) 

(220,335) 



Row-wise 

5. 

(197,437) (307,511) 

(252,501) 


120 

Row-wise 
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Fig. 3.9. a Some tiling arrangements for the rail road contact problem b,c Tiles used for unwrapping 
corresponding to Table 3.10-3.12 d Phase map used for unwrapping e Sample data points marked on the 
dark field image corresponding to Table 3.13. 
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Table 3.12 Unwrapping parameters and tiling plan 2 for phase unwrapping of rail road contact 
problem 


Tile 

Tile Coordinates 

Seed Point 

Seed Value 
(N) 

Threshold 

(tp) 

Seeding 

I. 

(40,94) (484,326) 

(159,236) 

2.00 

140 

Row-wise 

2. 

(87,324) (196,368) 

(188,334) 

3.00 

140 

Row-wise 

3. 

(195,324) (278,370) 

(240,359) 

3.00 

140 

Column-wise 

4. 

(276,324) (410,370) 

(327,368) 

3.00 

140 

Row-wise 

5. 

(172,366) (231,426) 

(221,371) 

4.00 

140 

Row-wise 

6. 

(229,369) (281,426) 

(245,374) 

3.00 

140 

Row-wise 

7. 

(280,368) (342,429) 

(293,398) 

3.00 

140 

Row-wise 

8. 

(205,423) (242,481) 

(232,472) 

6.00 

120 

Row-wise 

9. 

(240,424) (287,479) 

(252,477) 

3.00 

140 

Row-wise 

10. 

(284,427) (310,482) 

(289,475) 

4.00 

120 

Row-wise 

11. 

(194,480) (257,506) 

(239,501) 

4.00 

120 

Row-wise 

12 

(256,477) (321,505) 

(284,500) 

3.00 

120 

Row-wise j 


Table 3.13 Fringe orders at some salient points for the rail road contact problem unwrapped with 
different tiling arrangements shown in Tables 3.10-3.12 


No. 

Coordinates 

Unwrapped Fringe Order (N) 




1 

2 

1. 

(86,180) 

0.38 

0.38 

0.39 

2. 

(140,304) 

1.39 

1.40 

1.40 

3. 

(244,315) 

2.20 

2.21 

2.21 

4. 

(308,304) 

0.18 

0.19 

0.19 

5. 

(312,109) 

0.20 

0.21 

0.22 

6. 

(237,372) 

3.28 

3.29 

3.30 

7. 

(211,380) 

4.20 

4.21 

4.23 

8. 

(215,396) 

4.10 

4.10 

4.13 

9. 

(252,413) 

3.15 

3.16 

3.17 

IH 

(309,387) 

- 1.67 

2.33 

3.33 

11. 

(238,494) 

3.10 

3.12 

4.14 
































Fig. 3.10 Three dimensional fringe order plots for phase map of rail road contact problem (Fig, 3.9d) 
unwrapped with different tiling arrangements shown in Tables 3.10-3.12 respectively 
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Note that unwrapped fringe orders at some of the points are not matching 
back to phase map. These fringe orders are shown bold. Even with the arrangement 1 
these fringe orders are not improved. Only marginal improvement could be achieved by 
using arrangement 2. This is because phase map in these zones is not good and contains 
noise. Therefore, even with best possible tiling arrangement, if phase map contains 
noise and is not good, unwrapping will not be effective. 

3.9 Tiling for Problem of Spanner and Nut 



Fig. 3.1 1. Tiles for problem of spanner tightening the nut 

Tiling arrangement for the problem of spanner and nut is shown in Fig. 3.11. It is clear 
that as the model geometry looses symmetry and becomes more and more complex, the 
freedom in selection of tiles is reduced and more skills are required for the proper 
selection of tiles. One of the other problem in phase unwrapping comes into picture 
here. Disturbance in the image due to loading fixture can be seen. During unwrapping in 
and 5* tiles, this disturbance will propagate. Further, unwrapping procedure will 
completely be disturbed if the total fringe order at the seed points in neighboring tiles is 
supplied based on the results of 1®‘ or 5* tile. 

This requires a better selection of tiles for proper phase unwrapping. One of the 
other approach is shown in Fig. 3.12. Here tile 1 is divided in tiles 1 and 10. Also tile 5 
is divided in tiles 5 and 9. Now unwrap tiles 1 and 5 separately. And for all other tiles 
start unwrapping procedure from tile 6 and proceed in the sequence 7-8-10-4-9-2-3 or 
so. Care must be taken in selecting primary seed points for tiles 9 and 10 and should not 
be based on results of tile 1 or 5 as they will be unreliable due to noise in phase map. 



Fig. 3.12. Improved tiling arrangement for the problem of spanner tightening the nut 

Phase maps for the problem of spanner tightening the nut is shown in Fig. 3.13. 
Because of com.plexity in model geometry as well as loading condition, stress field is 
very complex. Also, since exact behaviour of isoclinic data is not known apriory, phase 
map generation is quite involved. In such problems it is very difficult to identify the 
exact interface of phase reversal. Figure 3.13a shows phase map generated with six-step 
phase shifting images using cos+sin expression. Since the phase map generated by 
cos+sin expression is not good, interactive approach was used for phase map 
generation. Phase map generated with interactive approach is shown in Fig. 3.13b. 



Fig. 3.13. Phase maps for Spanner problem: a Phase map obtained with cos+sin expression b 
Phase map obtained with interactive approach 

Unwrapping parameters used for unwrapping complete image are shown in 
Table 3.14. Table 3.15 shows the different tile coordinates and corresponding 
unwrapping parameters used to unwrap the phase map of Fig. 3.13b for the problem of 
spanner and nut with 512x512 image size. Corresponding tiles on a dark field image are 
shown in Fig. 3.14a. Table 3.16 shows the unwrapped fringe orders at some salient 




40 


points for each of unwrapping scheme. These points are shown, marked on a dark 
field image in Fig. 3.14b. 


Table 3.14 Unwrapping parameters used for unwrapping complete image for spanner and nut problem 


Tile 

Tile Coordinates 

Seed Point 

Seed Value 

(N) 

Threshold 

(tp) 

Seeding 

1. 

(0,0) (511,511) 

(98,158) 

1.00 

140 

Row-wise 


Table 3.15 Unwrapping parameters and tiling plan for phase unwrapping of spanner and nut problem 


Tile 

Tile Coordinates 

Seed Point 

Seed Value 
(N) 

Threshold 

(tp) 

Seeding 

1. 

(47,59) (129,262) 

(101,154) 

1.000 

140 


2. 

(88,6) (311,146) 

(96,119) 

0.674 

140 


3. 

(300,17) (498,148) 

(306,100) 

0.368 

140 

Row-wise 


(6,123) (73,315) 1 (64,160) j 0.494 i 140 ^Row-wise 

5. I (123,140) (336,199) (127,163) 0 88 140 Row-wise 


6. 

(333,146) (497,207) 

(361,148) 

0.400 

140 

Row-wise 

7. 

(237,195) (394,29) 

(295,198) 

1.152 

140 

Row-wise 

8. 

(388,200) (499,300) 

(432,207) 

0.498 


Row-wise 

9. 

(319,288) (49,380) 

(343,290) 



Row-wise 

10. 

(272,377) (495,446) 

(342,400) 

1.000 

140 

Row-wise 


11 . 


(195,428) (294,507) 


(248,489) 


1.000 


140 


Row-wise 


12 


(290,442) (441,505) 


(293,485) 


2.345 


140 


Column-wise 


Table 3.16 Fringe orders at some salient points for the spanner and nut problem obtained by unwrapping 
with and without tiling. 



(403,282) 0.15 


11. 

(426,279) 

0.21 

0.20 

12. 

(120,199) 

2.95 

2.94 
























































Fig. 3.14. Problem of Spanner and Nut a Tiles used for unwrapping corresponsing to Table 3.15. b 
Sample points of Table 3.16 marked on the dark field image 



Fig. 3.15 Three dimensional fringe order plots for phase map of Fig. 3.14b for the problem of spanner 
tightening a nut. a Unwrapped without tiling (SW isometric) b Unwrapped without tiling (NW isometric) 
c Unwrapped with tiles shown in Table 3.15 (SW isometric) d NW isometric view of(c) 
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f ToKu TIC U 1 ™^8e. Fringe orders unwrapped with scheme 

of Table 3.15 are shown in Fig. 3.15c and 3 15d i>cneme 

phase 4‘duiTotrg“ *■' 

order plot. This was effectively removed by the use Jf ming. 

3.10 Closure 


From the above drscussron it is clear that tiling procedure can be usefhl for a good 
number of praencal unvnapping problems. Also ttis can help in unwrapping in high 
ftmge graA^t zones as unwrapping threshold T, can be selected lower for particular 
Ule of htgh fringe gradient. Obviously the tiling concept calls for more user interaction 
Error propagation due to noise can be prevented by proper selection of tiles. 

understanding of fringe order variation from dark and bright field images 
and identification of phase reversal interfaces is important for generation of correct 
phase map, as no tiling can provide good results in phase unwrapping if, phase map 
Itself is not correct. Generated phase map should tally with corresponding dark and 

bright field images, if not, more refinement of phase map by interactive approach 
should be done [13]. 


Currently the software accepts only rectangular tiles. Integration of boundary 
identification routines with tile selection for accepting any generic tile shape can further 
reduce the calculational overheads. 



CHAPTER 4 


STRESS SEPARATION USING 
SHEAR DIFFERENCE TECHNIQUE 


4.1 Introduction 

Photoelasticity provides two independent equations for 2-D problem and five 
independent equations for 3“D problem. Since there are three unknown stress 
components for 2-D and six unknown stress components for 3-D problem, the equations 
obtained are not sufficient to determine the complete state of stress at a point. Hence 
with a view to obtain additional relations, several auxiliary methods are required for 
both 2-D & 3-D problems, which in conjunction with photoelastic analysis can provide 
the complete information for solving a general stress problem. 

Stress separation techniques, in a restricted sense, can be used to find the 
individual principal stresses and in a complete sense, they can provide all the 
components of the stress tensor. These methods can be either point-by-point or whole 
field. This is one of the most widely used methods in photoelasticity. However, it 
inherently bears the property of accumulation of errors. Tesar has modified the shear 
difference technique to minimize error accumulation [4]. The isoclinic parameter which 
is the chief source of error, is averaged along the three lines of grid and assumed to 
represent the isoclinic along the center line. 

Correct determination of isoclinic parameters in tiie field of photoelasticity is 
still a challenging task. Since, correct determination of isoclinic parameter is the basic 
requirement for success of shear difference technique, application of median filter 
during isoclinic determination is discussed first. 

4.2 Problems in Isoclinic Determination [7] 

The evaluation of isoclinic parameter is comparatively easier while determining it 
manually but the automated evaluation has several problems. The problems concerning 
the whole field isoclinic parameter evaluation can be roughly divided into the following 
parts. 
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1 . Ambiguity in experimentally evaluating the isoclinic parameter. 

2. Numerical inaccuracies. 

It is to be noted that calculations for isoclinic calculation involve inverse 
trigonometric operations and thus one has a multi-valued solutions, this poses a problem 
of correctly evaluating the isoclinic parameter. In most of the equations the values 
returned by the arctangent functions are in the range of-;c'4 to +7dA only. In a model, a 
priori one does not know whether the fast or slow axis lies the range -;z/4 to + 7 d 4 and 
whether the isoclinic determined, corresponds to oi direction or 02 direction. Thus the 
goal of isoclinic parameter evaluation should be such that it either gives cti or 
direction uniformly all over the domain. This problem could be solved [3] if theta is 
unwrapped in the range of - 71^2 to +n;/2. 

The problem of getting numerically correct value depends on the accurate 
recording and processing of the intensity information. In practical applications due to 
intensity digitization and quantization, only discrete values are used in the calculations. 
Also in the regions close to the isochromatic fringes, the error increases because of the 
ratios of small discrete values. Also there is a significant difference in the intensity 
maps of the images in theoretical and experimental conditions. These reasons lead to an 
error in the calculations. 

4.3 Median Filtering in Isoclinic Calculation 

Since shear difference is essentially a successive integration technique, a noise in either 
isochromatic or isoclinic parameter will disturb the integration procedure and successive 
results will loose their significance. Figures 4.1a and 4.1b show the isoclinic plot and 
theta field obtained by theoretically simulated images of disk under diametral 
compression for Brown and Sullivan algorithm. It can be seen that for all points where 
isochromatic fringe passes, the isoclinic determined has noise and oscillations. Causing 
the theta determined near an isochromatic to be unreliable and is called isochromatic 
interaction in determination of isoclinic. This suggests the noise removal in isoclinic 
determination. 

For noise removal in isoclinic, first the theta field is plotted as an intensity map 
of theta value ranging from —n/A to + 7 dA over the 0-255 gray levels. Then noise present 
in the theta field is removed by the use of median filter. Now, the new intensity values 
are read from the image buffer, converted back to theta value and put in the file 
containing isoclinic data {*.iso). Isoclinic plot and theta field obtained after application 
of above procedure on Fig. 4.1b are shown in Fig. 4.2a and 4.2b. 
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Fig. 4.1. a Isoclinic plot for the problem of disk under diametrical compression, obtained by 
application of Brown and Sullivan algorithm on theoretically simulated images b Corresponding theta 
field 



Fig. 4.2. Improved isoclinic data after application of median filter on Fig. 4.1b 

However, for experimental images, interaction of isoclinic and isochromatic 
does not cause noise but produces oscillation in theta near an full order isochromatic. 
Isoclinic plot and theta field obtained with experimental images are shown m Fig. 4.3a 
and 4.3b. Thus the above median filtering technique serve no purpose m isoclimc 
improvement for experimental images as median filter can only remove salt and pepper 
noise. As reported by Vinayak [7] isoclinic determined will be more reliable if images 
are taken at low loads, i.e. isochromatic interaction is reduced. Figure 4.3a and 43b 
shows the isoclinic plot and theta field obtained by low load experimental images of 
Brown and Sullivan algorithm. 
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Fig. 4.3. a Isociinic plot for the problem of disk under diametrical compression, obtained by 
application of Brown and Sullivan algorithm on experimental images b Corresponding theta field 



Fig. 4.4. a Isoclinic plot for the problem of disk under diametrical compression, obtained by 
application of Brown and Sullivan algorithm on Low Load experimental im.ages b Corresponding theta 


field 


It can be seen that theta determined with low load experimental images is better 
for the irmer regions of model. But, near the boundary, due to interaction with zeroth 
fringe order, isoclinic determined are bad. Thus, a judicious combination of above two 
results can provide a reasonably accurate isoclinic value for further calculations of stress 
fields by shear difference technique, which is discussed in next section. 
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4.4 Shear Difference Technique for Stress Separation 

4. 4.1 Conventional Method 

In the absence of body forces, the equation of equilibrium when applied to the 
problem of plane stress are, 

dx dy 
dy dx 



Fig. 4.5. Figure showing the grid needed for employing the shear difference method. The free 
boundary stress component is also shown. 

In the shear difference method, one of the equations in Eq. (4.1) is integrated to 
separate the stresses. Thus one can do stress separation for lines parallel to x-axis 
or y-axis. If in a problem, the stresses have to be separated along lines parallel to 
x-axis, the first of Eq. (4.1) is integrated to become. 


(gx) i~(<yx)^ 


‘ J Ph, 


' dy 


dx 


(4.2) 


Similarly for lines parallel toy-axis second of Eq. (4.1) leads 

I 

Using the finite difference approximation, Eq. (4.2) can be rewritten as, 


(Gy)j (Gy)i- J ^ 


(4.3) 




(4.4) 


where, (aj,, (o^), are the stress value at the start of the integration proced 
the starting point is chosen on a free boundary then (oi), is (if one has to do stress 
separation for lines parallel to x-axis) 
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(<Jx)i= gicos^6 = — — (4.5) 
h 

Once is calculated for a particular point on the grid, <Jy can be calculated as, 

(ay)^= cos 2B (4.6) 


On similar lines, integration of the equilibrium equations can be done for 
the };-direction. 

In above equations (o-x) u (cy ) , at the start of integration procedure chosen as 
boundary points where one of the component of stresses is zero. Thus isochromatic 
parameter at the boundary directly gives the value of the stress at the start of 
integration procedure. Xxy is then calculated at interior points and the summation 
proceeds in the stepwise maimer using fixed grid with each block measuring Ax by 
Ay. Axxy is obtained from the difference between the top and bottom values at each 
interior points. In above equations the convention for positive x is from left to right 
while for positive y is from bottom to top. Here in present work the direction of the 
scanning is taken from left to right and right boundary is assumed to be straight 
vertical line. Equations derived above are also based on left to right scanning i.e. 
by standard convention. 

The technique has come to be known as shear-difference technique because of 
the manner in which shear slope is calculated in Eq.{4.4) For performing the integration, 
one requires the data of both isochromatic and isoclinic fringe orders along the line of 
interest and for two more lines one above and the other below the line of interest. 

One of the basic problems associated with the technique is that the 
numerical errors accumulate as one proceeds along the line of interest. 
Nevertheless, the technique has been quite widely used because of its simplicity. 

4.4.2 Improvement by Tesar [4] 

One of the main sources of error accumulation in the conventional shear difference 
is the error in calculating the shear slope. It is well known that numerical 
differentiation is usually less accurate than numerical integration. Further, the 
calculation of the shear stress involves the use of isoclinic parameter w'hose 
measurement accuracy is generally poorer than that of the isochromatic fringe 
order. Tesar suggested that the shear stress be expressed in terms of the difference 
of principal stresses as, 


''xy 


(O’! -^2) 


sin 26 = 


^sinie 

2h 


2 


(4.7) 
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and be differentiated as a product of two functions. Substituting this into Eq. 
(4.2) gives, 

( CTx )j = ( 26 dx - \ai - a.jcos 26— dx (4.8) 

i dy 2 ■> dy 

Equation (4.8) can be approximated using the finite difference expressions as 

( <Jx )j 26 or CTi^cos 26— Ax. (4.9) 

/ 4y ^ i Ay 

The rest of the procedure is same as that of the conventional method. 

4.5 Results and Discussion 

Shear difference technique uses the isochromatic and isoclinic information for 
determination of normal stresses ( 0x and Oy) and shear stress (Xxy). Thus for all reasons, 
success of shear difference technique depends on correct determination of these 
parameters. 



Fig. 4.6. Pseudo fringe contours showing variation of stresses, determined using theoretical data by 
Tesar’s Modified Method a a, b Oy c Txy 
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Theoretical contours of Gx, Gy and Xxy for the problem of disk under diametral 
compression are shown in Fig. 4.6a, 4.6b and 4.6c respectively. These were obtained by 
first generating the isochromatic and isoclinic parameters theoretically, followed by the 
Tesar s modified stress separation technique on these data. Noise on the 
middle row is caused due to inaccurate generation of isoclinic data on that row. For 
stress separation, integration was carried out for each horizontal line from extreme left 
to extreme right. 



Fig. 4.7. Pseudo fringe contours showing variation of stresses, determined using data extracted from 
theoretically simulated images by Tesar’s Modified Method a Gx b Gy c 

Figure 4.7 shows the separated stresses obtained by carrying out the actual 
calculations on theoretically simulated images. It can be seen that, even for the 
theoretically simulated images, stress contours does not match to the actual one. This is 
because, calculation of theta even from simulated images is not accurate. One can see 
that, for each row in the image Gx and Gy plots are good until the first isochromatic 
fringe is reached. Since theta calculated at this point is not accurate, separated stresses 
contain error. This error is accumulated in subsequent integration and stress plots 
afterwards are not good. 

Comparison of above two case clearly indicates that error propagation due to 
numerical inaccuracies is negligible and main source of error in stress separation is 
inaccuracies in isoclinic calculations. This is further supported by the fact that, stress 
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contours obtained with data of simulated images show good match with theoreticallv 

calculated contours, till the fust isochromatic in a row is crossed, i.e. first inaccuracv in 
theta. 



Fig. 4.8. Pseudo fringe contours showing variation of stresses, determined using data extracted from 
experimental images by Tesar’s Modified Method a a* b Oy c i^y 

Separated stresses obtained with experimentally recorded images are shown in 
Fig. 4.8. Retardation calculation is carried out by six step phase shifting algorithm. 
Isoclinic parameters are calculated with Brown and Sullivan algorithm. These data were 
then put in a data file (*shd). Stress separation was carried out by Tesar’s technique. 

From the above two cases it can be seen that inaccuracies in calculation of 
isoclinic parameters does not affect the shear stress contours. This is because shear 
stresses are not determined by integration but is directly calculated for each point using 
isochromatic and isoclinic value at that point. Since at isochromatic fringes theta value 
become indeterminate, shear stress calculated at those points is also inaccurate, but they 
do not propagate. 


133624 
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4.6 Closure 

Problems in isoclimc determination and some techniques for isoclinic improvement is 
discussed. Then stress separation technique was presented. Since it was observed that 
t e c e source of error in stress separation is inaccurate determination of theta value, 

on 3.ccur3.tc C3.1cul3,tioii of th.ct 3 , is rccjuircd 



CHAPTER 5 


CONCLUSION AND SUGGESTIONS 
FOR FUTURE WORK 


5.1 Conclusion 

Unlike in conventional analysis, with advancement in digital photoelasticity data is 
available at pixel level. Thus the correct determination of fringe order values at each 
pixel in the domain is a crucial step. The main aim of this thesis is to develop an 
appropriate scheme, which will handle the phase map data efficiently for phase 
unwrapping, in the context of digital photoelasticity. The concept of tiling for effective 
phase unwrapping in the whole domain is presented. The applicability of the present 
scheme is verified for a variety of problems. The results of phase unwrapping are shown 
in the form of three dimensional fringe order plots. 

Accurate evaluation of isoclinic parameter 0 is still a challenging task in the 
field of digital photoelasticity. The evaluation of isocline parameter 6 from the plane 
polariscope arrangement gives reasonably good results. Also isoclinic parameter 
obtained at low loads in circular polariscope arrangement is comparatively better. So 
these methods are recommended for future use. However, further refinement m isoclinic 
parameter evaluation should be the focus of future research. 

4.2.Suggestion for Future Work 

1. Currently the software uses only rectangular tiles. This may be undesirable in 
certain class of problems due to complex geometry of the model. Use of boundary 
identification routines for generating tiles of any arbitrary shape can further improve 
the efficiency and effectiveness of phase unwrapping by tiling. 

2. The software in Windows environment made for photoelastic analysis stores 
unwrapped fringe order data and isoclinic parameter in the form of a float value 
written directly to a file in text format. This puts memory overheads in storing these 
data. The use of binary file format to store these data can considerably reduce the 
memory required for storage of these data. 

3. From results of in-plane shear stress txy, it can be seen that the isoclinic parameter 
does not have much effect on it. This tells that if some smoothening of G is done 
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before using it for stress separation, there is a higher probability of satisfactory 
results. 

Since the isoclinic parameter Q, plays an important role in the end results, there is a need 
for improvement of it. 
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